library(texreg)
load('nes.congress.RData')

#Constructing Binned Levels of Congruence
nes.congress$cov2c.3bins <- cut(nes.congress$cov2c, breaks=c(0, .25, .63, 1), labels=FALSE)
nes.congress$cov2c_92.3bins <- cut(nes.congress$cov2c_92, breaks=c(0, .25, .52, 1), labels=FALSE)

#Subsetting to Races with Competition
nes.congress <- subset(nes.congress, nes.congress$runopposed==0)

#Function to fit regression with cluster-robust standard errors
fit.w.robust <- function(model,cluster.var){
	
	robust.se <- function(model, cluster){
		require(sandwich)
		require(lmtest)
		M <- length(unique(cluster))
		N <- length(cluster)
		K <- model$rank
		dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
		uj <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
		rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
		rcse.se <- coeftest(model, rcse.cov)
		return(list(rcse.cov, rcse.se))
	}
	
	out <- model$model
	clustervar <- mapply(paste,cluster.var,out[,grep(pattern=cluster.var, x=names(out))],sep="")
	vcov <- robust.se(model, clustervar)
	model$se <- vcov[[1]]
	model$coeftest <- vcov[[2]]
	return(model)
}

##############
#Within District
within.district <- lm(I(congress.vote=="Dem") ~ cov2c + I(PID=="Dem") + cov2c*I(PID=="Dem") + factor(distyear), data=nes.congress)
within.district <- fit.w.robust(model=within.district,cluster.var="distyear")

#Within District + Controls
within.district.controls <- lm(I(congress.vote=="Dem") ~ cov2c + I(PID=="Dem") + cov2c*I(PID=="Dem") + lpop + ldensity + ld2 + ld3 + ld4 + urban2 + curb2 + curb3 + curb4 + curb5 + school2 + school3 + linc + s20l + s65m + black + fema + factor(age) + factor(gender) + factor(edu) + factor(urban) + factor(inc) + factor(distyear), data=nes.congress)
within.district.controls <- fit.w.robust(model=within.district.controls,cluster.var="distyear")
coeftest(within.district.controls, vcov=within.district.controls$se)

#SE on main coefficient
within.district.coefficient <- within.district.controls$coefficients[length(within.district.controls$coefficients)]
within.district.se <- sqrt(within.district.controls$se[length(within.district.controls$coefficients),length(within.district.controls$coefficients)])
within.district.upper <- within.district.coefficient+2*within.district.se
within.district.lower <- within.district.coefficient-2*within.district.se
within.district.coefficient; within.district.lower; within.district.upper

#Within District + Controls + binned congruence
within.district.controls.binned <- lm(I(congress.vote=="Dem") ~ factor(cov2c.3bins) + I(PID=="Dem") + factor(cov2c.3bins)*I(PID=="Dem") + lpop + ldensity + ld2 + ld3 + ld4 + urban2 + curb2 + curb3 + curb4 + curb5 + school2 + school3 + linc + s20l + s65m + black + fema + factor(age) + factor(gender) + factor(edu) + factor(urban) + factor(inc) + factor(distyear), data=nes.congress)
within.district.controls.binned <- fit.w.robust(model=within.district.controls.binned,cluster.var="distyear")
binned.coef.table <- coeftest(within.district.controls.binned, vcov=within.district.controls.binned$se)

marginal.binned.coef <- binned.coef.table[dim(binned.coef.table)[1], 1]
se.binned.coef <- binned.coef.table[dim(binned.coef.table)[1], 2]
marginal.binned.coef; marginal.binned.coef-2*se.binned.coef; marginal.binned.coef+2*se.binned.coef

##############
#Redistricting
redistricting <- lm(I(congress.vote=="Dem") ~ cov2c_92 + I(PID=="Dem") + cov2c_92*I(PID=="Dem") + factor(countyid) + factor(stateyear), data=nes.congress)
redistricting <- fit.w.robust(model=redistricting,cluster.var="countyid")

redistricting.model.controls <- lm(I(congress.vote=="Dem") ~ cov2c_92 + I(PID=="Dem") + cov2c_92*I(PID=="Dem") + factor(age) + factor(gender) + factor(edu) + factor(urban) + factor(inc) + I(incumbent.party=="Dem") + age_rep + freshman + tenure + died + retired + powercom + status + incrun + factor(countyid) + factor(stateyear) + I(incumbent.party=="Dem")*age_rep + I(incumbent.party=="Dem")*freshman + I(incumbent.party=="Dem")*tenure + I(incumbent.party=="Dem")*died + I(incumbent.party=="Dem")*retired + I(incumbent.party=="Dem")*powercom + I(incumbent.party=="Dem")*status + I(incumbent.party=="Dem")*incrun, data=nes.congress)
redistricting.model.controls <- fit.w.robust(model=redistricting.model.controls,cluster.var="countyid")

#Redistricting + Controls + binned congruence
redistricting.model.controls.binned <- lm(I(congress.vote=="Dem") ~ factor(cov2c_92.3bins) + I(PID=="Dem") + factor(cov2c_92.3bins)*I(PID=="Dem") + factor(age) + factor(gender) + factor(edu) + factor(urban) + factor(inc) + I(incumbent.party=="Dem") + age_rep + freshman + tenure + died + retired + powercom + status + incrun + factor(countyid) + factor(stateyear) + I(incumbent.party=="Dem")*age_rep + I(incumbent.party=="Dem")*freshman + I(incumbent.party=="Dem")*tenure + I(incumbent.party=="Dem")*died + I(incumbent.party=="Dem")*retired + I(incumbent.party=="Dem")*powercom + I(incumbent.party=="Dem")*status + I(incumbent.party=="Dem")*incrun, data=nes.congress)
redistricting.model.controls.binned <- fit.w.robust(model=redistricting.model.controls.binned,cluster.var="countyid")

#############################################
#Produce Tables / Graphs for in-text material
congruence.sequence <- seq(from=0,to=1,by=.05)

#Marginal Effects for Continuous Congruence
interaction.number <- grep(pattern='cov2c',x=names(within.district.controls$coefficients))[2]
me.plot <- within.district.controls$coefficients[3] + within.district.controls$coefficients[interaction.number]*congruence.sequence
me.standard.error <- sqrt(within.district.controls$se[3,3] + congruence.sequence^2*within.district.controls$se[interaction.number,interaction.number] + 2*congruence.sequence*within.district.controls$se[3,interaction.number])

interaction.number.redistrict <- grep(pattern='cov2c',x=names(redistricting.model.controls$coefficients))[2]
me.plot.redistrict <- redistricting.model.controls$coefficients[3] + redistricting.model.controls$coefficients[interaction.number.redistrict]*congruence.sequence
me.standard.error.redistrict <- sqrt(redistricting.model.controls$se[3,3] + congruence.sequence^2*redistricting.model.controls$se[dim(redistricting.model.controls$se)[1]-8,dim(redistricting.model.controls$se)[1]-8] + 2*congruence.sequence*redistricting.model.controls$se[3,dim(redistricting.model.controls$se)[1]-8])

#Point Estimates for Binned Congruence
within.district.controls.binned.coefs <- c(within.district.controls.binned$coefficients[4], within.district.controls.binned$coefficients[4] + within.district.controls.binned$coefficients[grep(pattern=':I\\(PID', x=names(within.district.controls.binned$coefficients))])

redistrict.controls.binned.coefs <- c(redistricting.model.controls.binned$coefficients[4], redistricting.model.controls.binned$coefficients[4] + redistricting.model.controls.binned$coefficients[grep(pattern=':I\\(PID', x=names(redistricting.model.controls.binned$coefficients))])

#Combined Marginal Effects Plot
pdf(file="votebycongruence-combinedplot.pdf", height=4, width=10)
par(mar=c(4.1,6,2,2), mfrow=c(1,2))
plot(y=me.plot, x=congruence.sequence, las=1, xlab='Congruence', type='n', ylim=c(.325,.7), ylab='Increase in Pr(Vote) For Co-Partisan', main='Within-District')
polygon(x=c(congruence.sequence, rev(congruence.sequence)), y=c(me.plot+2*me.standard.error,rev(me.plot-2*me.standard.error)), border=NA, col='gray75')
points(y=within.district.controls.binned.coefs, x=c(.125,.44,.815), cex=.8, pch=16,col='black')
lines(x=congruence.sequence, y=me.plot, lwd=2.5)

plot(y=me.plot.redistrict, x=congruence.sequence, las=1, xlab='Congruence', type='n', ylim=c(.325,.7), ylab='Increase in Pr(Vote) For Co-Partisan', main='Redistricting')
polygon(x=c(congruence.sequence, rev(congruence.sequence)), y=c(me.plot.redistrict+2*me.standard.error.redistrict,rev(me.plot.redistrict-2*me.standard.error.redistrict)), border=NA, col='gray75')
points(y=redistrict.controls.binned.coefs, x=c(.125,.385,.76), cex=.8, pch=16,col='black')
lines(x=congruence.sequence, y=me.plot.redistrict, lwd=2.5)
dev.off()


#Typical Change In Congruence
residualized.congruence <- lm(cov2c ~ lpop + ldensity + ld2 + ld3 + ld4 + urban2 + curb2 + curb3 + curb4 + curb5 + school2 + school3 + linc + s20l + s65m + black + fema + factor(age) + factor(gender) + factor(edu) + factor(urban) + factor(inc) + factor(distyear), data=nes.congress)$residuals
typical.congruence.change <- sd(residualized.congruence)

residualized.congruence.redistrict <- lm(cov2c_92 ~ factor(age) + factor(gender) + factor(edu) + factor(urban) + factor(inc) + I(incumbent.party=="Dem") + age_rep + freshman + tenure + died + retired + powercom + status + incrun + factor(countyid) + factor(stateyear) + I(incumbent.party=="Dem")*age_rep + I(incumbent.party=="Dem")*freshman + I(incumbent.party=="Dem")*tenure + I(incumbent.party=="Dem")*died + I(incumbent.party=="Dem")*retired + I(incumbent.party=="Dem")*powercom + I(incumbent.party=="Dem")*status + I(incumbent.party=="Dem")*incrun, data=nes.congress)$residuals
typical.congruence.change.redistrict <- sd(residualized.congruence.redistrict)

#Parametric Bootstrap for counterfactual from model
set.seed(seed=1130)
par.boot <- rnorm(n=2000,mean=within.district.controls$coefficients[interaction.number],sd=sqrt(within.district.controls$se[interaction.number,interaction.number]))
par.boot <- par.boot*typical.congruence.change
typical.change <- typical.congruence.change*within.district.controls$coefficients[interaction.number]
typical.change
quantile(x=par.boot, probs=c(.025, .975))

par.boot.redistrict <- rnorm(n=2000,mean=redistricting.model.controls$coefficients[interaction.number.redistrict],sd=sqrt(redistricting.model.controls$se[length(redistricting.model.controls$se[3,]),length(redistricting.model.controls$se[3,])]))
par.boot.redistrict <- par.boot.redistrict*typical.congruence.change.redistrict
typical.change.redistrict <- typical.congruence.change.redistrict*redistricting.model.controls$coefficients[interaction.number.redistrict]
typical.change.redistrict
quantile(x=par.boot.redistrict, probs=c(.025, .975))

##################
#Tables of Results
##################
names(redistricting.model.controls$coefficients)[2] <- 'cov2c'
names(redistricting.model.controls$coefficients)[length(names(redistricting.model.controls$coefficients))-8] <- 'cov2c:I(PID == "Dem")TRUE'

within.district.controls.se <- sqrt(diag(within.district.controls$se))
within.district.controls.pvalue <- coeftest(within.district.controls, vcov=within.district.controls$se)[,4]

redistricting.controls.se <- sqrt(diag(redistricting.model.controls$se))
names(redistricting.controls.se)[2] <- 'cov2c'
names(redistricting.controls.se)[length(names(redistricting.controls.se))-8] <- 'cov2c:I(PID == "Dem")TRUE'

redistricting.controls.pvalue <- coeftest(redistricting.model.controls, vcov=redistricting.model.controls$se)[,4]
names(redistricting.controls.pvalue)[2] <- 'cov2c'
names(redistricting.controls.pvalue)[length(names(redistricting.controls.pvalue))-8] <- 'cov2c:I(PID == "Dem")TRUE'

texreg(list(within.district.controls,redistricting.model.controls), omit.coef=c('factor|lpop|ldensity|ld|urban|curb|school|linc|s20l|s65m|black|fema|(Intercept)|age|freshman|tenure|majparty|died|retired|powercom|status|party|rankmem|chair|close|incrun|runopposed'), stars=c(0.05), caption="Pr(Vote for Politician) by Level of Congruence", custom.model.names=c("Within-District","Redistricting"), override.se=list(c(sqrt(diag(within.district.controls$se))),c(sqrt(diag(redistricting.model.controls$se)))), override.pval=list(c(coeftest(within.district.controls, vcov=within.district.controls$se)[,4]), c(coeftest(redistricting.model.controls, vcov=redistricting.model.controls$se)[,4])))
#Note: function incorrectly reports a warning message, the correct cluster-robust SE's are reported

###################
#Robustness Checks for District Competitiveness
###################

redistricting.model.controls.competition.cook <- lm(I(congress.vote=="Dem") ~ cov2c_92 + I(PID=="Dem") + cov2c_92*I(PID=="Dem") + factor(age) + factor(gender) + factor(edu) + factor(urban) + factor(inc) + I(incumbent.party=="Dem") + age_rep + freshman + tenure + died + retired + powercom + status + incrun + factor(countyid) + factor(stateyear) + I(incumbent.party=="Dem")*age_rep + I(incumbent.party=="Dem")*freshman + I(incumbent.party=="Dem")*tenure + I(incumbent.party=="Dem")*died + I(incumbent.party=="Dem")*retired + I(incumbent.party=="Dem")*powercom + I(incumbent.party=="Dem")*status + I(incumbent.party=="Dem")*incrun + factor(competition.cook), data=nes.congress)
redistricting.model.controls.competition.cook <- fit.w.robust(model=redistricting.model.controls.competition.cook,cluster.var="countyid")

names(redistricting.model.controls.competition.cook$coefficients)[2] <- 'cov2c'
names(redistricting.model.controls.competition.cook$coefficients)[length(names(redistricting.model.controls.competition.cook$coefficients))-8] <- 'cov2c:I(PID == "Dem")TRUE'

redistricting.model.controls.competition.cq <- lm(I(congress.vote=="Dem") ~ cov2c_92 + I(PID=="Dem") + cov2c_92*I(PID=="Dem") + factor(age) + factor(gender) + factor(edu) + factor(urban) + factor(inc) + I(incumbent.party=="Dem") + age_rep + freshman + tenure + died + retired + powercom + status + incrun + factor(countyid) + factor(stateyear) + I(incumbent.party=="Dem")*age_rep + I(incumbent.party=="Dem")*freshman + I(incumbent.party=="Dem")*tenure + I(incumbent.party=="Dem")*died + I(incumbent.party=="Dem")*retired + I(incumbent.party=="Dem")*powercom + I(incumbent.party=="Dem")*status + I(incumbent.party=="Dem")*incrun + factor(competition.cq), data=nes.congress)
redistricting.model.controls.competition.cq <- fit.w.robust(model=redistricting.model.controls.competition.cq,cluster.var="countyid")

names(redistricting.model.controls.competition.cq$coefficients)[2] <- 'cov2c'
names(redistricting.model.controls.competition.cq$coefficients)[length(names(redistricting.model.controls.competition.cq$coefficients))-8] <- 'cov2c:I(PID == "Dem")TRUE'

texreg(list(redistricting.model.controls,redistricting.model.controls.competition.cook,redistricting.model.controls.competition.cq), omit.coef=c('factor|lpop|ldensity|ld|urban|curb|school|linc|s20l|s65m|black|fema|(Intercept)|age|freshman|tenure|majparty|died|retired|powercom|status|party|rankmem|chair|close|incrun|runopposed'), stars=c(0.05), caption="Pr(Vote for Politician) by Level of Congruence", custom.model.names=c("Redistricting","Redistricting+Cook Competition","Redistricting+CQ Competition"), override.se=list(c(sqrt(diag(redistricting.model.controls$se))),c(sqrt(diag(redistricting.model.controls.competition.cook$se))),c(sqrt(diag(redistricting.model.controls.competition.cq$se)))), override.pval=list(c(coeftest(redistricting.model.controls, vcov=redistricting.model.controls$se)[,4]), c(coeftest(redistricting.model.controls.competition.cook, vcov=redistricting.model.controls.competition.cook$se)[,4]),c(coeftest(redistricting.model.controls.competition.cq, vcov=redistricting.model.controls.competition.cq$se)[,4])))
#Note: function incorrectly reports a warning message, the accurate cluster-robust SE's are reported